{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from k_means import KMeans\n",
    "from rv import RandomVariable\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "class MultivariateGaussianMixture(RandomVariable):\n",
    "    \"\"\"\n",
    "    p(x|mu, L, pi(coef))\n",
    "    = sum_k pi_k N(x|mu_k, L_k^-1)\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self,\n",
    "                 n_components,\n",
    "                 mu=None,\n",
    "                 cov=None,\n",
    "                 tau=None,\n",
    "                 coef=None):\n",
    "        \"\"\"\n",
    "        construct mixture of Gaussians\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        n_components : int\n",
    "            number of gaussian component\n",
    "        mu : (n_components, ndim) np.ndarray\n",
    "            mean parameter of each gaussian component\n",
    "        cov : (n_components, ndim, ndim) np.ndarray\n",
    "            variance parameter of each gaussian component\n",
    "        tau : (n_components, ndim, ndim) np.ndarray\n",
    "            precision parameter of each gaussian component\n",
    "        coef : (n_components,) np.ndarray\n",
    "            mixing coefficients\n",
    "        \"\"\"\n",
    "        super().__init__()\n",
    "        assert isinstance(n_components, int)\n",
    "        self.n_components = n_components\n",
    "        self.mu = mu\n",
    "        if cov is not None and tau is not None:\n",
    "            raise ValueError(\"Cannot assign both cov and tau at a time\")\n",
    "        elif cov is not None:\n",
    "            self.cov = cov\n",
    "        elif tau is not None:\n",
    "            self.tau = tau\n",
    "        else:\n",
    "            self.cov = None\n",
    "            self.tau = None\n",
    "        self.coef = coef\n",
    "\n",
    "    @property\n",
    "    def mu(self):\n",
    "        return self.parameter[\"mu\"]\n",
    "\n",
    "    @mu.setter\n",
    "    def mu(self, mu):\n",
    "        if isinstance(mu, np.ndarray):\n",
    "            assert mu.ndim == 2\n",
    "            assert np.size(mu, 0) == self.n_components\n",
    "            self.ndim = np.size(mu, 1)\n",
    "            self.parameter[\"mu\"] = mu\n",
    "        elif mu is None:\n",
    "            self.parameter[\"mu\"] = None\n",
    "        else:\n",
    "            raise TypeError(\"mu must be either np.ndarray or None\")\n",
    "\n",
    "    @property\n",
    "    def cov(self):\n",
    "        return self.parameter[\"cov\"]\n",
    "\n",
    "    @cov.setter\n",
    "    def cov(self, cov):\n",
    "        if isinstance(cov, np.ndarray):\n",
    "            assert cov.shape == (self.n_components, self.ndim, self.ndim)\n",
    "            self._tau = np.linalg.inv(cov)\n",
    "            self.parameter[\"cov\"] = cov\n",
    "        elif cov is None:\n",
    "            self.parameter[\"cov\"] = None\n",
    "            self._tau = None\n",
    "        else:\n",
    "            raise TypeError(\"cov must be either np.ndarray or None\")\n",
    "\n",
    "    @property\n",
    "    def tau(self):\n",
    "        return self._tau\n",
    "\n",
    "    @tau.setter\n",
    "    def tau(self, tau):\n",
    "        if isinstance(tau, np.ndarray):\n",
    "            assert tau.shape == (self.n_components, self.ndim, self.ndim)\n",
    "            self.parameter[\"cov\"] = np.linalg.inv(tau)\n",
    "            self._tau = tau\n",
    "        elif tau is None:\n",
    "            self.parameter[\"cov\"] = None\n",
    "            self._tau = None\n",
    "        else:\n",
    "            raise TypeError(\"tau must be either np.ndarray or None\")\n",
    "\n",
    "    @property\n",
    "    def coef(self):\n",
    "        return self.parameter[\"coef\"]\n",
    "\n",
    "    @coef.setter\n",
    "    def coef(self, coef):\n",
    "        if isinstance(coef, np.ndarray):\n",
    "            assert coef.ndim == 1\n",
    "            if np.isnan(coef).any():\n",
    "                self.parameter[\"coef\"] = np.ones(self.n_components) / self.n_components\n",
    "            elif not np.allclose(coef.sum(), 1):\n",
    "                raise ValueError(f\"sum of coef must be equal to 1 {coef}\")\n",
    "            self.parameter[\"coef\"] = coef\n",
    "        elif coef is None:\n",
    "            self.parameter[\"coef\"] = None\n",
    "        else:\n",
    "            raise TypeError(\"coef must be either np.ndarray or None\")\n",
    "\n",
    "    @property\n",
    "    def shape(self):\n",
    "        if hasattr(self.mu, \"shape\"):\n",
    "            return self.mu.shape[1:]\n",
    "        else:\n",
    "            return None\n",
    "\n",
    "    def _gauss(self, X):\n",
    "        d = X[:, None, :] - self.mu\n",
    "        D_sq = np.sum(np.einsum('nki,kij->nkj', d, self.cov) * d, -1)\n",
    "        return (\n",
    "            np.exp(-0.5 * D_sq)\n",
    "            / np.sqrt(\n",
    "                np.linalg.det(self.cov) * (2 * np.pi) ** self.ndim\n",
    "            )\n",
    "        )\n",
    "\n",
    "    def _fit(self, X):\n",
    "        cov = np.cov(X.T)\n",
    "        kmeans = KMeans(self.n_components)\n",
    "        kmeans.fit(X)\n",
    "        self.mu = kmeans.centers\n",
    "        self.cov = np.array([cov for _ in range(self.n_components)])\n",
    "        self.coef = np.ones(self.n_components) / self.n_components\n",
    "        params = np.hstack(\n",
    "            (self.mu.ravel(),\n",
    "             self.cov.ravel(),\n",
    "             self.coef.ravel())\n",
    "        )\n",
    "        while True:\n",
    "            stats = self._expectation(X)\n",
    "            self._maximization(X, stats)\n",
    "            new_params = np.hstack(\n",
    "                (self.mu.ravel(),\n",
    "                 self.cov.ravel(),\n",
    "                 self.coef.ravel())\n",
    "            )\n",
    "            if np.allclose(params, new_params):\n",
    "                break\n",
    "            else:\n",
    "                params = new_params\n",
    "\n",
    "    def _expectation(self, X):\n",
    "        resps = self.coef * self._gauss(X)\n",
    "        resps /= resps.sum(axis=-1, keepdims=True)\n",
    "        return resps\n",
    "\n",
    "    def _maximization(self, X, resps):\n",
    "        Nk = np.sum(resps, axis=0)\n",
    "        self.coef = Nk / len(X)\n",
    "        self.mu = (X.T @ resps / Nk).T\n",
    "        d = X[:, None, :] - self.mu\n",
    "        self.cov = np.einsum(\n",
    "            'nki,nkj->kij', d, d * resps[:, :, None]) / Nk[:, None, None]\n",
    "\n",
    "    def joint_proba(self, X):\n",
    "        \"\"\"\n",
    "        calculate joint probability p(X, Z)\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        X : (sample_size, n_features) ndarray\n",
    "            input data\n",
    "\n",
    "        Returns\n",
    "        -------\n",
    "        joint_prob : (sample_size, n_components) ndarray\n",
    "            joint probability of input and component\n",
    "        \"\"\"\n",
    "        return self.coef * self._gauss(X)\n",
    "\n",
    "    def _pdf(self, X):\n",
    "        joint_prob = self.coef * self._gauss(X)\n",
    "        return np.sum(joint_prob, axis=-1)\n",
    "\n",
    "    def classify(self, X):\n",
    "        \"\"\"\n",
    "        classify input\n",
    "        max_z p(z|x, theta)\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        X : (sample_size, ndim) ndarray\n",
    "            input\n",
    "\n",
    "        Returns\n",
    "        -------\n",
    "        output : (sample_size,) ndarray\n",
    "            corresponding cluster index\n",
    "        \"\"\"\n",
    "        return np.argmax(self.classify_proba(X), axis=1)\n",
    "\n",
    "    def classify_proba(self, X):\n",
    "        \"\"\"\n",
    "        posterior probability of cluster\n",
    "        p(z|x,theta)\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        X : (sample_size, ndim) ndarray\n",
    "            input\n",
    "\n",
    "        Returns\n",
    "        -------\n",
    "        output : (sample_size, n_components) ndarray\n",
    "            posterior probability of cluster\n",
    "        \"\"\"\n",
    "        return self._expectation(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x1 = np.random.normal(size=(100, 2))\n",
    "x1 += np.array([-5, -5])\n",
    "x2 = np.random.normal(size=(100, 2))\n",
    "x2 += np.array([5, -5])\n",
    "x3 = np.random.normal(size=(100, 2))\n",
    "x3 += np.array([0, 5])\n",
    "x_train = np.vstack((x1, x2, x3))\n",
    "\n",
    "x0, x1 = np.meshgrid(np.linspace(-10, 10, 100), np.linspace(-10, 10, 100))\n",
    "x = np.array([x0, x1]).reshape(2, -1).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARMAAAD8CAYAAABUzEBbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXmcVMW5v5+afaZnBlxAcQGJEhDEBRA1UeKGATeMQUVz\n3a77z2sSE5O4xCWauOVqYlyuYsQl7hsBFUVxCSGCiIoLi4KALCKbLMIMDDP9/f1R3TM906dneobu\nnp6Z9+lPffqcOlV16pyhv1S99VaVk4RhGMa2ktPaFTAMo31gYmIYRkowMTEMIyWYmBiGkRJMTAzD\nSAkmJoZhpISUiIlzboxzbqVz7rOYuO2dc2845+ZFvrdLkHeYc+5z59x859yVqaiPYRiZJ1Utk0eA\nYQ3irgTelNQLeDNyXg/nXC5wLzAc6Auc7pzrm6I6GYaRQVIiJpImA982iB4BPBo5fhQ4KSDrYGC+\npAWSqoCnI/kMw2hj5KWx7J0kLY8cfwPsFJBmV2BJzPlS4KCgwpxzFwIXAoRCoYF9+vRJYVUNw4jl\ngw8+WC2pS3PypFNMapEk59w2+e1LGg2MBhg0aJBmzJiRkroZhhGPc+6r5uZJ52jOCudcN4DI98qA\nNMuA3WPOd4vEGYbRxkinmIwHzo4cnw2MC0jzPtDLOdfTOVcAjIrkMwyjjZGqoeGngKlAb+fcUufc\necCtwFDn3Dzg6Mg5zrldnHMTACRVA/8DTATmAM9KmpWKOhmGkVlSYjORdHqCS0cFpP0aODbmfAIw\nIRX1MAyj9TAPWMMwUoKJiWEYKcHExDCMlGBiYhhGSjAxMQwjJZiYGIaREkxMDMNICSYmhmGkBBMT\nwzBSgomJYRgpwcTEMIyUYGJiGEZKMDExDCMlmJgYhpESTEwMw0gJJiaGYaQEExMjc3z5JZx9Nuy1\nFxx9NLz1VmvXyEghGVmd3jD44gs48EDYtAlqarywTJ0KDzwA//VfrV07IwWktWXinOvtnJsZEzY4\n537ZIM3hzrn1MWmuS2edjFbi2mth40YvJFEqKuCXv6wfZ7RZ0toykfQ5sD/UbgW6DBgbkPTfko5P\nZ12MVubf/4ZwOD6+shKWLYPu3TNfJyOlZNJmchTwpaRmb+5jtAN23jk4vqYGtgvc095oY2RSTEYB\nTyW49gPn3CfOuVedc/0yWCcjU1x1FZSU1I8rKoJTT4Wystapk5FSMiImkQ22TgSeC7j8IdBd0r7A\n3cA/E5RxoXNuhnNuxqpVq9JXWSM9nHIK3HgjlJZ68SgqghEjvAHWaBc4aZu2AE7uJs6NAC6VdEwS\naRcBgyStTpTG9hrOYsJhP3JTUhJsB6mshAULfLdnhx0yXz8jKZxzH0ga1Jw8mermnE6CLo5zbmfn\nnIscD47UaU2G6mWkkkmTYNddYdAg6N3bf3/VwERWXAz9+pmQtEPS7mfinAsBQ4GLYuIuBpB0PzAS\nuMQ5Vw1UAqOUieaSkVoWLvTdloqKuriPPoLDD/c+JTkx/29VV8PYsTBunBeVCy6AffZp+h4SvP02\nPPIIVFXBz34Gxx1Xv2yj9ZDU5sLAgQNlZBm/+52Uny/5n3xdKCuTHntMevVVafFiqapKGjJECoX8\n9dxcqbhYeuSRpu9xxRV1+UAqLZVGjZLC4fQ/XwcDmKFm/i4zYjNJNWYzyUJOPx2efjo+PjfXtxxK\nSmDLFhgwAGbOrN+CAX99xQpvoAVvV1m2DPr3h86dvR1mv/1g8+b6+UIhmDABhgxJz3N1ULLZZmK0\nRySYPx8WL/ZzbUKh+DQ1NbB1K6xf74Vg2rR4IQHIy4MpU2DtWt812mcfOOEE6NYNrr8eJk4MrsOm\nTfDKKyl9LKNlmJh0ZKqq4M9/hu9/H/bYA377W/+jT4apU6FnT99a6N0b7rgDCgubtl8EecGCF6ZQ\nyNtBpk71oz5RAbrjDpg717dyGpKfD+XlydXZSCvWzenIHHecN2hWVvrzwkIvEB9/DAUFifOtXAl7\n7unn2gThnA9R60Yy5OXBBx/A4MG+O5SozIZiVFwMs2d7MTRSRku6OTZruKPywQfwzjt1QgL+R7x0\nKbzwgreBRJkzB2691YvMgAGw445+RCYRURHJyYkXE+ega1dvH2mY54ILvKgEiUmsMIVCvpVSXQ0P\nPWRCkiWYmHRUpk8PbjVs3OhtF1ExmTrV20O2bPH2j88+84LQmJhEibYioq2UggJvaC0rixeTmhpv\nmA2FvB0kETk53n/l8svhyCPNFT+LMJtJR2X33X0roCHFxfC97/njcNivNVJRUbdMQE2NF5Lm+HZI\nvjXzi1/A++/HO7JFycuDP/whfg5PLOGwF6IRI0xIsgwTk47KsGHQqVO8KOTlwVlneQEYNcoP0QYR\nDgcbRBPx7bfw3Xfw2muJ7Sjl5XDppd6OM2KEb9E0JCfHDxcbWYeJSUclL8+vMXLggd7wWlTkR3Xe\nfBO6dPFdnQkTEufffnvo1Sv5+4XDcP/9cN99ibtIZ5wBd9/tF1JascLXrbi4fpqiIrjmmuTva2QM\ns5l0ZPbYw/t9rFzpfUF22aWuNfDaa8H+IOCHYy+7zIvOb3+bOF0Qc+f61kXDUZmcHPjrX+vHFxd7\nG0penm/V9O/vxWa//Zr1mEZmMDEx/OhKQzp18gbToJGVrl1962DLFu8DsmyZ91lJhqhHbMP0iVZh\ncw7uvNOP9NgcnKzG/jodna1b/SS9hj4jZ5yR+Me7di08/rh3fZ8xw4+s7LVXcjaU6mpvYI2O6hQV\nBdtGolRUwMsvm5C0Aewv1JG57z7fVenf339fcEFdi2G33eDvfw/OV1HhWyTgWzAHHOCXFWhMFGJZ\nt86LSUWFF7PGRCg3N/GSj0ZWYd2cjsqLL8JvflPf3vHEE/7He//9/nzIEN9yaDi5DmDVKj/57pe/\nhMmT431DunTxaYLYcUdYsyY579iCAj/CY2Q91jLpqPzpT/GG08pKePTRuvhddkm82PPq1bDvvt5Q\n21BIrrjCG3VfecWPFMVy/vnwzTfwr38l5ydSVeV9T779NrnnMloNE5OOytKlwfHO1f1wc3L8Gq0l\nJXU2i2hXJhz2BtiGrYvf/MZPHgQ49lj45z/rBOX88+HBB33r57DD4NVX6wtKcbHvXsV2e2pqvCgd\nc0zy83yMVsHEpKNy8MHBNo6SEj/tP8oJJ/huzKmnertIY/aN3/wGbr/dC80NN/jWybBhXlAuvdQL\nCcDNN3sv2B/+sL6gVFZ6kWu4KdfWrX5I+YMPtumRjTTT3NWUsiHYSmsp4NNP/UplztWtXFZSIo0Z\nkzjPq69KnTrFr6YG0qWX+jQ1NdK55/q4fv2kFSvql3H55f5az57SokU+bsoUv+JaULnRUF4uPfdc\net6FEQctWGnNWiYdlX328Q5rP/mJt40cfDA89xyce27iPAccEOx3AnX7BU+fDo895o9nzfKT8Vau\n9Oe/+hX85S/+eOFC7z8CvoXS1Mzfqiq/WPUjj/hh6XXrknlKI5M0V32aG4BFwKfATALUDnDA34D5\nwCfAgKbKtJZJK3L55fXXYY2GwYOldet8mmeeqd/S6NpVGjCgfvqRI6WtW336m25qvFVSXCwdcIBU\nVORbU6WlPu7FF1vvPbRzgn6rTYVMicmOjVw/Fng1IioHA+81VaaJSYb58kvpP/+RNmzwizc/+qhU\nWNg8QWmpkIB03HFeSIJEZvXq1nsv7Zi2KiYPAKfHnH8OdGusTBOTDLFmjXTYYf5H26mT/77lFn+t\nT5/gH/7gwdLGjT7NlVfGX99jjzohueWWpoUEpLw8HxrGh0LS3//eOu+mndMSMcmEzUTAJOfcB865\nCwOu7wosiTlfGomrh20P2gqceiq8917deqyVlfDHP3rbym67BefZd9+6haUXLoy/vm5d3dDz/vvH\n+6EkImimcXR42sgOmqs+zQ3ArpHvrsDHwJAG118GDo05fxO/Pai1TFqTZcuCuxbR/WoKCuLjzz+/\nLn901CYo9O1bN8rz6qvBXabYEDviFBuKiupGhIyUQja2TCQti3yvBMYCgxskWQbsHnO+WyTOaE3W\nrPFLDQSxaVP8rN+zzqrzI4kdtQG/dcXPf17nozJ7NhxxRH0/lKjPS5DviwKc1QoL/bonPXo067GM\n9JFWMXHOhZxzZdFj4BjgswbJxgNnOc/BwHpJy9NZLyMJvve9+otNxxL0447On5k61a9LEuWUU+CN\nN+Cuu+DJJ+sLyh/+4I+HDfOzjnNy4OyzveOcc4lnCjvn16U99VQ47TTYaSfvUPfII+Yl25o0tynT\nnAB8D9+1+RiYBVwTib8YuDhy7IB7gS/xQ8iNdnFk3ZzM8L//G9yVSRSOOEKqqPB5H3jAd01OOaXO\n2Lp5s/+OjvIcc4xUWenj7rmnrtuy++7SCy9Iv/61dNpp3pEu6H6DBnmjcE5OXVxJiXTNNZl/V+0Q\nsnE0Jx3BxCQDJBqtaSycf36doEycWCckN95Yf9j4zTfrhOTee+uXkZMj9e/v8++wQ/B9Cgu9mATt\nbVxUJK1fn/n31c4wMTFSx557Nl9MQDryyDpBkbyQxA4bRwVFiheSWEFIZPzNy/NOcHvvHXy9vFya\nNi3z76ud0RIxMXd6I5jTT09+2DaWt96C44/3E/OuuQauu67u2vTpfvbvrFlw222J1ymJbqfREOdg\n6FD46CPo0yfYWFtV5bfxMDKObQ9qBPPdd37OTNCSjskQCnlBaIkfSH6+nykcVObdd/v5Q9On+xGh\n2DVZioq8WI0b1/x7GvWw7UGN1FFW5qf8jxvnt8To1Mn/mG+4IXjltYZs2hS8yVcydO0Ky5fHLzId\nDsOhh/rjwYP96NAll/g1aSU4+WQYPbpl9zS2GWuZdFS++w7mzfNdgi5dks+3ww7Jr3pWUJD8qvWx\n5OXFd3NCIT8UPGZM/fhw2K/cFhU7IyW0pGViNpOOhgS//733zTjiCC8mo0Yl19oA7weS7MLRLRES\niBcS5/wyk0ELXOfk+CUUTEhaHROTjsbf/+69UysrYcMGb9MYP97vA5wMN9wAe++d2Ds2luZ2cxI5\nqYVCfhFq2+4iq7G/Tkfj9tuDF5J+7LHkjKXl5fDxx/DMM9C9e+OtlOrq5FsxXbv6LTeCloXMyUm8\nsLWRNZiYdDRWrw6OD4eTH7XJy/MrtH38sZ/52xgNbXKJWhcrV3oX+4brv4Ifoh46NLm6Ga2GiUlH\n49BDg1sLO+3kNyNvDp07++UIioqSzxO0DWiU6HBwTo5vAZWV+aUaJ01KrltltComJh2N227z23pG\nuxPO+Yl1992XfJcklj339H4fsQbQnJyWlRXLk0/CxImweLFfI8XIeszPpKPRty98+KHfbmLaNOjV\nC66+Gg46qOVl3nuvXzj6/vu9f8mBB/oh3IabcyVLXp53PrPWSJvC/EyM1CN5p7JPP22ZB+zPfuZX\noDdaDfMzMbID5/wcnR/9qPl5Cwt9l8toc5iYGOmhrAzmzGlenqIiv2BSeXl66mSkFbOZGOmhpgaW\nLEl8PTfXt2Byc31rpKoKTjwRfv3rzNXRSCkmJkZ6yM31o0aJfFfOP98bbmfO9DOT99/fL91otFms\nm2Okj1GjguPz8+u8XQcOhJEjTUjaAeleUHp359zbzrnZzrlZzrm4CSDOucOdc+udczMj4bqgsow2\nyN13e9tJQ4qL4YwzMl8fI62ku2VSDfxaUl/81p+XOuf6BqT7t6T9I+HGNNfJyBRFRX5NlP3398sR\nFBZC797w9ts216YdklabifyWFcsjx9855+bgd+ubnc77GllEr15+mcXly71RNtFOgEabJ2M2E+fc\nHsABwHsBl3/gnPvEOfeqc65fgvy2PWhbpls3E5J2TkbExDlXCrwA/FLShgaXPwS6S9oXuBv4Z1AZ\nkkZLGiRpUJfmrAxmGEZGSLuYOOfy8ULyhKQXG16XtEHSxsjxBCDfObdjuutlGEZqSfdojgMeAuZI\nujNBmp0j6XDODY7UaU0662UYRupJt9PaD4EzgU+dczMjcVcD3QEk3Q+MBC5xzlUDlcAotcXZh4bR\nwUn3aM4U/F7CjaW5B7gnnfUwDCP9mAesYRgpwcTEMIyUYGJiGEZKMDExDCMlmJgYhpESTEwMw0gJ\nJiaGYaQEExPDMFKCiYlhGCnBxMQwjJRgYmIYRkowMTEMIyWYmBiGkRJMTAzDSAkmJoZhpAQTE8Mw\nUoKJiWEYKcHExDCMlJCJ1emHOec+d87Nd85dGXDdOef+Frn+iXNuQLrrZBhG6kn36vS5wL3AcKAv\ncHrA9qDDgV6RcCHwf+msk2EY6SHdLZPBwHxJCyRVAU8DIxqkGQE8Js80oLNzrlua62UYRopJt5js\nCiyJOV8aiWtuGtse1DCynDZjgLXtQQ0ju0m3mCwDdo853y0S19w0hmFkOekWk/eBXs65ns65AmAU\nML5BmvHAWZFRnYOB9ZKWp7lehmGkmHTv6FftnPsfYCKQC4yRNMs5d3Hk+v3ABOBYYD5QAZybzjoZ\nhpEe0r3XMJIm4AUjNu7+mGMBl6a7HoZhpJc2Y4A1DCO7MTExDCMlmJgYhpESTEwMw0gJJiaGYaQE\nExPDMFKCiYlhGCnBxMQwjJRgYmIYRkowMTEMIyWYmBiGkRJMTAzDSAkmJoZhpAQTE8MwUoKJiWEY\nKcHExDCMlGBiYhhGSjAxMQwjJaRt2Ubn3J+BE4Aq4EvgXEnrAtItAr4DaoBqSYPSVSfDMNJHOlsm\nbwD7SNoX+AK4qpG0R0ja34TEMNouaRMTSa9Lqo6cTsPvh2MYbZ5Fi+Cyy2DwYDj3XJg1q7VrlB2k\nfXX6CP8NPJPgmoBJzrka4AFJo4MSOecuxG9sTvfu3dNSScNoitmz4eCDYfNm2LoVPvwQnnsOJkyA\nIUNau3atyzaJiXNuErBzwKVrJI2LpLkGqAaeSFDMoZKWOee6Am845+ZKmtwwUURkRgMMGjRI21Lv\nbKaaat7lXWqo4Qf8gEIKW7tKRgxXXAEbN4Ii/wJramDTJrj4Yi80HZltEhNJRzd23Tl3DnA8cFRk\nf5ygMpZFvlc658YCg4E4MWlrCDGf+eSTzx7skVSeKUzhJE5iK1tr457iKY7l2DTV0oj+q3QuufST\nJ9flieWLL6CiAkpKUle3tkbabCbOuWHAb4ETJVUkSBNyzpVFj4FjgM/SVadMMY1p9KQn+7M/felL\nP/oxl7mN5tnABo7lWNawhg0xn5GM5Gu+zlDNOw6rV8MZZ0BRERQUwIknwtKlTefr3Dk4Pj/fl9OR\nSedozj1AGb7rMtM5dz+Ac24X51x0h7+dgCnOuY+B6cArkl5LY53SzipWMZShfMVXVFBBJZXMYQ5D\nGMIsZnEN13ARFzGe8dRQU5tvLGMJE44rL0yYJ3kyk4/Q7qmpgUMPheefh6oqqK72No/Bg6GysvG8\nQ4fGxxUUwNlnQ16mLJBZStoeX9JeCeK/xu8tjKQFwH7pqkM62MxmXuEVVrGKIQyhL33rXX+cx+uJ\nBPguz0Y2MoABCLGVrTzJkxzIgUxkIvnks5a1VFNNQ7awhTWsSeszdTQmToSvv/YG1Cg1NfDdd/Ds\ns14Ygvj6a3j66fj4mhq4/vr01LUt0cG1tGlWspKHeZh5zKMHPfgrf2UrW2t/+KdyKmMYQ06kkbeE\nJVQS/99bw7iNbGQ603mCJziHcziKo2rLiCVEiMM4jDu5k/GMZytbWc96qqjiOI7jKq6iK13T8OTt\nl7lzYcuW+PiNGxsf5n366WB7SU4O/Nd/wVlnwamnQnFx6urappDU5sLAgQOVCWZohspUpiIVCSEn\nJxp8QgrpCT1Rm2esxqpUpXHpEn2Gamht3vN0nkIK1St7uIZrL+2lYhXH5c1Xvrqpm9ZoTUbeR3vh\n1VelsjLJS0NdKC2VHn00cb7f/15yLj5fbP4ePaRvvmm6Ds88I+2zj7TddtLQodIHH6Ts8VICMEPN\n/F22ujC0JGRKTPbW3kkJwhANqc2zVVs1UANrBQihQhUqRzmBeUtVqkt0iRZqocIKa6zG6ngdr6Ea\nqpEaqXKVN3rvPOXpRt2YkffRXqiulvbeW8rLqy8GO+wgVVQkzvevf0kFBYnFBHyZ55zT+P3/9jep\npKR+vpIS6aOPUvuc20JLxMQm+iVgBStYwIKk0m6hrs2cRx6Tmcy1XEsf+rAjO1JDTaBxFXx35wEe\nYD/2YzKTWcxijuAI1rOe8YxnAxsavXc11TzMw8k/mEFuLrz0UvzoS0UF/N//wYYN8PLLMGlSfbvK\nhAnePtIY1dUwdmzd+YYN8Nln/ht8edde6+8VS2Wlj2/TNFd9siFkomXyrb5VgQqabJWUqER36+7A\nMh7TY/W6LU19cpSjIhUldd+G3Z1qVaf9nbQnrr9eKiyMb1kUFEhFRVJ5uQ/bby9NnSqtW+fjG2uV\nRENxsXTffdLFF/s8ZWX++xe/kBYujG+VRMPOO7f2W6mDFrRMzACbgO3YjkM4hClMiRudySWXGmoo\npZR92ZcLuCCwjDu4g01sSvqeYcJsZnOz6+pwfMu3dKFLs/N2VCZMCDbCVlX5780xf4Zhw3xLpalW\nSZTKSvj5z30rJbasBx+E8vJgIy7AHnskV362YmLSCE/wBIdxGKtZXSsoB3EQP+JHtf4kx3M8ueTW\n5lke+fSmN9/ybUbqWUghnUngTWXUsnYt3HsvvPaaH+ZNlnAYXn21fpenKarjR/mpqID77oPzzoMx\nY+p3dUpK4Lrrki8/GzExaYRd2ZV5zOMt3mIRixjEIA7ggMC0G9nI6ZzOa7xGmDBC7MIu5JEX6D+S\nKgoo4AquIJ/8tN2jPbB6NRxwgP/enKDx51xwqyEc9naPVLB2LfzlL97BbfRo39opL4c77oDhw1Nz\nj1ajuf2ibAiZGs1pDifr5MARmxzlqFCFtcf5yo+ziTi5hKM9DT+xw9N5ytPNullhhVv78bOe3/0u\n2EYC3jZSXCz17h1sz8jNla6+OnH+5oT99qur0+bN0sqVUk1N672XRGA2k9ZhLWt5iZcCR2yEGMYw\nvuALCihgJCM5hmN4lEdZzGK6051VrOI//IcVrMDhqKEGEdyxPpmT6U9/DuVQjuCIQEc3I56XXw62\nkZSW+q7PD37gbRbDhvlRnNgWSk0N3H13YltHc7jyyrrjwkLo0o7MXCYmLeQLvuAu7mI2s+lLXxzB\n006FmMIUtrCFTWziS77kLu7ip/yUWcxiAn6aUpgwueSSQw496MGXfBlY3lu8xfM8n7bnao8E2S+i\nbN0Khx0GPXr48xNOgDffjBeO777zQ8pFRYm7SU1RWAgHHdSyvG0BE5MWMJnJDGc4VVTVrj9SRVXC\n9OtYV2vA3Rj5jGZ0XOujJvJprKy1rGUQg5jAhDg3eiHe4R1e4iXKKONMzmQvAqdIdRhWrvStjmXL\ngq9XV9cZQquq4Ne/9jaSIGpqkh/RCaJnTx/aLc3tF2VDaE2bSVhhfV/fT8q+kc5PJ3XSVbpKy7Vc\nklSjGp2qU2v9WvKVr2IV6zE91mrvKhs45ZR4T9fY4Jz0gx/4tFdc0bi9o3t3qVev4Gvdukl9+ybO\nm5cnffpp676L5oC506eXsML6Wl8rX/nbLAZ91EeTNVnX6Jq4az3UQ2/qTd2qWwPnA8UaYDurs+Zr\nvsZrfKCDXLGKtU7rWuV9tTbhcNPu71EDa2Wl1KVL4jQ/+pG0aZNUVSWdeGL9awMHSmvXekPqmWcm\nLuPII6UNG1r7rSSHiUmaCCusu3SXdtAOjf64myMky8PLa8u/TbfVE5IFWlB7bbRGN3nPwRqsM3RG\n4LVylesFvVBb3mzN1p/0J92iW/SFvsjoe8w04XDjrZJoyM/3c3JychoXkiixghIVkihNCcrw4T5d\nZaV0883S97/vWzs33dT4vKBMY2KSJu7W3c1yi09WSGYsm6Gq6ipJXlBihWTm8pnaVOX/BScjKCM1\nMjBNgQp0uS5XhSp0k25SsYqVp7zabtBf9JeMvstMc9JJvuXRmJCcdJI0d64/bkxIHnpI+vOf/XFV\nlR9ujgrJ88/7WcVS04LSr58fio6dgVxcLB18cPYME5uYpIGwwtpROwb+gJ2cOqlT0nNpYoXk9fmv\nq/iPxRrx1IhaQVmv9ZKkqUumqvyWch356JFJC8oFukAlKgm8FlJIXdW13kzm6KdIRVqkRRl7n5lm\n2TJp99398gDgWx/O+R9vWZmfPbxypZ8z03B5gf796wtJ9Prtt9e/x/PP17WArrrKx9XUSMcem1hQ\ngkJpqfTaaxl9PQlpiZiYk0ITbGVrwpXO8sjjaZ7mMz6r51KfiFu4hZ3dzsxbM48RT4+gsrqScZ+P\n45TnTmFrzVbKKWfa0mn8+PEfs2HLBt5a+Bb/Pe6/AbiACziEQxKWvT3b80t+SRFFcb4nm9jEKlYF\nzvsRYhzjmqx7W2WXXWDePO9t+vvfw+OPwzvvwF//6mf3vvGG90a97LL44eD+/esWiH7yybrrv/0t\n3H67P37uORg1qm74+Z//9MPIOTl+GcjmsHEjTJsWH798uV8Vf8AAOPnk4DRZQXPVJ9kA3AAsA2ZG\nwrEJ0g0DPgfmA1cmU3amuzm7aJfA//H30T6SpCVaErh4UcPPKTpFNWHfjr1y0pXiBmrDEY8coT9N\n/pPKbymvjev6566atXKWJOkTfZKwq1WiEk3RFEnSAi1QrnKTaCfVtUz+pr9l9H1mA//5j7TXXn42\nb35+sL3EOenBB336TZukI46of71bt/rne+8tLY+Ywl5/PflZxtEQCvkRpZtu8vddt05avNivsxI1\nJDvnvXSfeSa974ds6uZExOSKJtLk4vch/h5QAHwM9G2q7EyLyT/0j7guRLGK9bJeluQXROqszkn9\neM/SWQkFJTY0FJIu6hJYXkghnaEz6rnUR933kxWTJVqS0ffZ2nz1VV23p6nQlKCkSkjAd5VCIX/P\nUMi7+Y8YEWxE3nFHv8hTumiLYnIIMDHm/CrgqqbKbo2h4ef0nHqrt4pUpP7qr1f0Sr3r9+iehDaL\nxgTlvHHnxQlJ3o15+nSFd0oIEpJCFaqP+ug0naZH9Iju1J26VbdqtmZLks7QGXHD1/nK1wE6QEUq\nUqEKVRShVrt9AAAR/klEQVT5PKAHMvsis4Crr05uyDhWUJ5/3uf95pvgNP/+t7/+/vstE5IddgjO\nl8h4HApJ8+en7x21REzSbTO5zDn3iXNujHNuu4DruwJLYs6XRuLicM5d6Jyb4ZybsWrVqnTUNRAh\nXuIlRjOaTnTiRm5kKlPjNsa6lEt5mIfpRz860YkjOZLruI4C4jdT+Rf/YqVbCcDAXQbGXS8vLKfX\n9r0AeId3WEXd8+aRxwpWMIc5HMdxXMzFXM3VXMu1DGQgV3M1d3EX3+N7lFFGAQWUUcae7MkbvMFs\nZnMrt3Ibt/EFX3Ch33G1QzFvXt26JckwYAAcdZQ/fv314DQTJ/rv/faDY45Jrtzo3JwnnoDddw92\n05fi48DbaLYL+kW1Js1Vn9gATMJvmtUwjMDviZOL35vnT8CYgPwjgb/HnJ8J3NPUfTPZMrlaV9ez\nVRSrWP3UTxWq7xQQVlgf6SNdqAvVQz20p/bUzbpZD+vhejaM2OHf6KhNUDcndpQn1g8lRzl6T+9p\nlVYFjs6UqETTNV3VqtYETdAdukMTNMFWYovh7rsTr3bWMMT6kcSO2gSF6ChPkGNbUDjhhLqh4P33\nD05TWOhHnmLjCgp83nRCNnVz6t0E9gA+C4jP6m7Oci0PtD+EFNKDerA23Wf6TL3UK27otljFOlJH\n6hk9IyfXqJAU/bFIfe7pk5SgFKpQl+myQINsjnL0K/2q0ed6US+qr/qqRCUaoAGaqInpe4lZyIYN\n0m67Nf1j79UrsZD07i09/LD0wx8mFpTDD09cdk6OX9YxStAi0+BXu7/hhrqlJIuKpKOP9sbZdJJV\nYgJ0izm+HHg6IE0esADoSZ0Btl9TZWdKTF7QCwlXhx+hEZKkClVoB+0QmCYqPAfpICH0vHzHe96a\neQlHbRoaZc8ae1ZtfQ7WwbXl5ikv4ejOSI1M+EyP6/FAY/JryhIHhwzx4INNe8eefXZd+kMOqYuP\nNbY2NMr26FHnMv/HP9YJR8OyS0qk6dPryq+q8iIRCvn0oZDUqZM0Y4a/vm6dH4FatCgTbyf7xOQf\nwKfAJ8D4qLgAuwATYtIdC3yBH9W5JpmyMyUmkzVZZSqL+7HmKleX6BJJ0lN6KjBNbCsiOtJzrs6t\nLfuily6KE5IoUUHpdEsnTV/q/8XN13x1Uqd6ZecpL+E9x2lc3POEFdau2jUwz77aN41vMvs477zG\nhQS82Dz7rE+/dq104IH1heTrr/13VFB69PDOb5L07ru+JbHddr4FU1bmQ2mpb100dHyTvPv/O+94\nEXroodadx5NVYpLOkCkxqVGNeqhH3CpoJSrRJ/pEknSn7mzUA7ZMZTpAB9Se/1w/ry3/d2/8rt7w\n7y/0i9pRnj+884d6QrK7do8rO5EwINRTPeOeZ7M2J/SiLVBBul9nVnHddcHu80GCMnasz7N2bd0G\nW6+/7lsXscPGixf746iQ7L679N13Pq6iQnrhBemxx+rEKJsxMUkD8zVffdRHIYVUrnKVqUxP6ana\n6+/pvUbn7eyoHTVFU+q1ImIFRao//Bs7bBy9f5CQoOAdBmOvPa/n9ZHqdnaqUU2jm4G1daqq/OS5\n7t29H8a559a1HhqyaFHyw8OxLRSpvh9JrB+KVCck+fl1ItQWMTFJE2GF9Zk+07t6V5u1Oe7a8Tq+\n0S7HhbpQ3dStXnxUUIL8SKKC0piQJPMpV7lCCmmgBmqVVulrfZ2wnuUqz+g7TQc/+Un9kY+8PO+l\nmshYecMNyYlJtKw77pAeeCDeH8Q56cYb/dai5eU+7pBD/Hm6DaXpwsSklVimZQl/pI19eqpnwi5S\nd3VPykU/mU++8jVcw7VBGxLer4/6tPZr3CbmzIkfQo0aOv/61+A8K1YkXnZgW0NRkbePhELSm29m\n9l2kgpaIiU30SwEf8mGgc1pTLGRhwiUaF7OYSioDrw1mMIUUJn2frWzlTd6khhpO4IS4vCWU8Ct+\nlXzFs5CPPvIT9hpSUQH/+U9wnq5d4eyz01OfzZv9xL1Nm+C44/zGXI3x0Udw5JEQCkH37qlbwDqT\nmJhsI2/wBqdxGhVUNJ24CQopbHL2cYgQYxnLOZxDMcX1rhVTTDnlgflyyGEjGxnDGA7lUIopphOd\nKKKIC7mQ8zl/m+vfmuyxR/DarYWF0Lt34nzN8YRtKZs3w0UXJb4+Z45f1Prtt734LVniV7GPXcm+\nTdDcpkw2hGzp5qzTupQtmpSjHB2mw/S4Hm90lOYKXaF1WqfN2qwzdIaKVKQylalYxbpe1+sCXRDY\n5eqhHvUmA87XfL2lt7RSK1vxDaaOcFjad9/4EZrSUmlJgnmM4XBq9sJJttuzeXNwPX72s+DuVlGR\ntH59+t5ZY2A2k8zymB5TqUoT/vDzla8iFSV0fGv4KVCB9tW+ulW3Bk4adHIKKaRSlep23a6X9JK+\n0lf6WB/rO/kxyK/1tXbWzrX2ljzlqUQlmqRJrfy20s+qVdLxx/tRmoICv8DztGmJ04fDja/ClspQ\nUpLY4ax37+A85eXSRx8F50k3LRET2+piG9jEprhNzaP8iB9xEiexmtXcwR1JlVdFFbOZTTXVHM/x\njGc8QO2iRkK1G6H/jt9RSilCPMuz7Mu+AHSjG7OZzWhG8zZv04teXMZlfJ/vb+vjZj077ggvveTt\nFFu2wPbbN57eORg61E/ea9hFysmJj3POx7d0u4uddgqO79MHvvgi3kZSVeXtJ22G5qpPNoRsaZks\n0ILAyXYhhfSO3tF6rW/RiEyBCnSOztHn+lx91KfJ9CUq0Wqtbu3X0SZZsEDafvs61/qcHKlzZ++F\nGvVYLSyUjjnGj8q0pFtUUiJdc03iOkyfHj8vp7hYOueczL2HhmCjOZmlJz25gisooaR2R78QIU7k\nRIYwhA/5sEUbildRxXM8xxrWsAM7JJXHdvlrGQWRQbhoqyAc9i2CLl38Bl5TpsDChX6JgW++qUvf\nnPKvvRZuvLF+fHU1rFvn73vggfDii7Dnnn7XwJISb7B94IFtf76M0lz1yYaQLS2TKFM0RefpPJ2p\nM/WKXqk1dH6qT5NeMCnoc7yO1026qUkjb77ydbsCJnsYTfL//l/whL+ysniD6fvve7+RZFskeXnS\nXXfVL6OmxrdSSku9XWennbyLfZRNm9K7glqyYAbY7KO/+geuyVqsYvVRH+2tvROKRIEKVKQidVM3\nhRRq1H1+uIarUpWt/bhtjj32CBaCsjLpk0/qpw2HvWdrMnN6wLv0r1vnBeTFF6Wf/lTq0ye+q1RS\nIr30Uus8fyJMTLKIsMK6STepRCW1IlCgAuUqV9tpOx2uwzVd07VGa5q0qxSrWJfrcl2v63WwDg7c\nUbBYxbpIF7X2Y7cpXn45sQdsYWHwhLy1a4M9bfPy/IZaxcW+xXHiiX4GcTjstyhtqkUzYEDGH79R\nTEyyiDt1Z2AXJ9pKyVGOSlSiZ/SMJmhC7ZBvIkGJLhEQVjhunk/0U6QibdGWVn7ytsHmzXXzaBqG\n3Ny6nfcasnBh4lXaunePT//OO8mt6rb99ml93GbTEjExA2yauIVbAr1io0PJYcJUUMFFXMTRHM0K\nVnA7tyd0k4+W5XAJvW3DhNnIxhQ9Qfvm3XcTX+vcGZ56KvhaaWnioeFOneLjxo71Xq1Nsd9+TafJ\ndkxM0oAQq1mdVNoqqpjDHEKEuJiL6UKXuDSFFHIap9WeH8qhtaNHsXSjG9uRbasMZydB83iiHHJI\nsDCA92U59FDIbzBIV1ICP/95fPqpU5uuS0kJ3Hxz0+myHROTNOBwSTuJbWELnelcm+8f/IMQodqJ\ngyFC7MEe/Ibf1Oa5jdsopbR2Ho/DUUIJ93JvoMgY8RxySPAwbygE553XeN4nn4R+/Xza8nIoKvIT\nBoPyzZ2buJzSUvjhD/2uggcf3Lz6ZyPOd4/aFoMGDdKMGTNauxqN8iqvMpKRTU4A3IEd4loxS1jC\nQzzEQhZyJEdyGqdRRFG9NPOZz83czHu8Ry96cTVXM5hm7kfZwXn3XRg+3PuWVFd7D9czz4T77/fH\njSHBzJmwdKnfCmPXwA1avBfu2rXx8bm53lO3MPnJ3xnFOfeBpEHNypMuMXHOPQNE52t2BtZJ2j8g\n3SLgO6AGqE7mAdqCmAC8zdtcy7XMYhYb2ECY+v7ZDseN3Mjv+X0r1dDYuNHvD/ztt3D00dC3b2rL\nv+QSGDOm/uzk3Fw44gjfIslWskpM6t3EuTuA9ZJuDLi2CBgkKTkjA21HTGI5kROZxKTaNUpyyaUr\nXZnL3ITLBhhtn3XrfFdm8WJviC0p8V2jqVOze95NS8Qk7RP9nHMOOBU4Mt33ag1mMpNJTKKcckYy\nku0Jnl32PM9zK7cymtFUUskJnMDN3GxC0s7p3Bk++QQmTPDfe+0FJ52Uvd2bbSHtLRPn3BDgzkQq\n55xbCKzHd3MekDS6qTKzoWUixIVcyJM8STXV5JGHwzGOcRzFUa1aN8PYVjLeMnHOTQJ2Drh0jaRx\nkePTgQSj9gAcKmmZc64r8IZzbq6kyQH3uhD8xrjds6B9+Aqv8BRP1RpYo8sv/pSfspKVLVrG0TDa\nMtskJpKObuy6cy4POBmI3527roxlke+VzrmxwGAgTkwiLZbR4Fsm21DtlPAwD9euLRKLEJOZzNE0\n+moMo92Rbj+To4G5kpYGXXTOhZxzZdFj4Bj8xudZT6JFkYC4URvD6AikW0xG0aCL45zbxTk3IXK6\nEzDFOfcxMB14RdJraa5TSjiLswgRiosXYghDWqFGhtG6pHU0R9I5AXFf4/cXRtICoE3OSjiJkziB\nE3iJl6iggkIKcTie4qk4BzPD6AjYGrAtJIccnuRJpjGNiUykM50ZxSh2DrRHG0b7x8RkG3A4Dol8\nDKOjYxP9DMNICSYmhmGkBBMTwzBSgomJYRgpwcTEMIyUYGJiGEZKMDExDCMlmJgYhpESTEwMw0gJ\nJiaGYaQEExPDMFKCiYlhGCnBxMQwjJRgYmIYRkowMTEMIyWYmBiGkRJMTAzDSAnbJCbOuVOcc7Oc\nc2Hn3KAG165yzs13zn3unPtxgvzbO+fecM7Ni3xvty31MQyj9djWlsln+H1x6u1z45zri1+Zvh8w\nDLjPOZcbkP9K4E1JvYA3I+eGYbRBtklMJM2R9HnApRHA05K2SFoIzMdvrhWU7tHI8aPASdtSH8Mw\nWo90LSi9KzAt5nxpJK4hO0laHjn+Br+PTiCx24MCW5xzbWKzrmayI7C6tSuRJtrrs7XX5+rd3AxN\nikmS+wlvM5LknEu47Wfs9qDOuRnN3VS5LdBenwva77O15+dqbp4mxaSp/YQTsAzYPeZ8t0hcQ1Y4\n57pJWu6c6wasbMG9DMPIAtI1NDweGOWcK3TO9QR64bf/DEp3duT4bCBlLR3DMDLLtg4N/8Q5txQ4\nBHjFOTcRQNIs4FlgNvAacKmkmkiev8cMI98KDHXOzcNvcn5rkrcevS31zmLa63NB+302e64ITkpo\npjAMw0ga84A1DCMlmJgYhpES2oyYbKvrflvBOXeDc26Zc25mJBzb2nXaFpxzwyJ/l/nOuXbl4eyc\nW+Sc+zTyd2r2UGq24Jwb45xbGeu71ZKpLm1GTNh21/22xF8k7R8JE1q7Mi0l8ne4FxgO9AVOj/y9\n2hNHRP5ObdnX5BH8byeWZk91aTNikgLXfSPzDAbmS1ogqQp4Gv/3MrIISZOBbxtEN3uqS5sRk0bY\nFVgSc57Idb8tcZlz7pNI87Mtz6Ruj3+bWARMcs59EJnu0Z5IeqpLlHTNzWkRmXLdb20ae07g/4Cb\n8P9QbwLuAP47c7UzmsGhkpY557oCbzjn5kb+l29XNDXVJUpWiUmaXfezhmSf0zn3IPBymquTTtrc\n36Y5SFoW+V7pnBuL79a1FzFp9lSX9tDNSdZ1v00Q+cNF+Qne8NxWeR/o5Zzr6ZwrwBvKx7dynVKC\ncy7knCuLHgPH0Lb/Vg1p9lSXrGqZNIZz7ifA3UAXvOv+TEk/ljTLORd13a8mxnW/jXK7c25/fDdn\nEXBR61an5Uiqds79DzARyAXGRKZatAd2AsY658D/jp6U9FrrVqllOOeeAg4HdoxMj7keP7XlWefc\necBXwKlNlmPu9IZhpIL20M0xDCMLMDExDCMlmJgYhpESTEwMw0gJJiaGYaQEExPDMFKCiYlhGCnh\n/wOO19+3quxqfQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f7a3497e9b0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "gmm = MultivariateGaussianMixture(n_components=3)\n",
    "gmm.fit(x_train)\n",
    "p = gmm.classify_proba(x_train)\n",
    "\n",
    "plt.scatter(x_train[:, 0], x_train[:, 1], c=p)\n",
    "plt.scatter(gmm.mu[:, 0], gmm.mu[:, 1], s=200, marker='X', lw=2, c=['red', 'green', 'blue'], edgecolor=\"white\")\n",
    "plt.xlim(-10, 10)\n",
    "plt.ylim(-10, 10)\n",
    "plt.gca().set_aspect(\"equal\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# training data\n",
    "x1 = np.random.normal(size=(100, 2))\n",
    "x1 += np.array([-5, -5])\n",
    "x2 = np.random.normal(size=(100, 2))\n",
    "x2 += np.array([5, -5])\n",
    "x3 = np.random.normal(size=(100, 2))\n",
    "x3 += np.array([0, 5])\n",
    "x_train = np.vstack((x1, x2, x3))\n",
    "\n",
    "x0, x1 = np.meshgrid(np.linspace(-10, 10, 100), np.linspace(-10, 10, 100))\n",
    "x = np.array([x0, x1]).reshape(2, -1).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARMAAAD8CAYAAABUzEBbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXecVNX5/99ne6VJERWsBCwgymoUjQiCBbuxxoKaxJ4v\nUWMUSwL6NRpbNPaC+jUhYAX9IQhYEjEqClgAC0WxIIKiIOzC1s/vjzOzO7NzZ3dnd2br857Xee0t\n55459w73wznPec5znCQMwzCaSlpLV8AwjPaBiYlhGEnBxMQwjKRgYmIYRlIwMTEMIymYmBiGkRSS\nIibOuUedc2udc4sjjnVzzs1xzi0L/e0a59rDnXOfOueWO+euSkZ9DMNofpLVMnkcOLzWsauAVyT1\nA14J7UfhnEsH7gWOAHYDTnPO7ZakOhmG0YwkRUwkvQ78UOvwscD/hbb/Dzgu4NJ9geWSPpNUBkwJ\nXWcYRhsjI4Vl95K0OrT9LdArIM+2wFcR+18DPw8qzDl3HnAeQH5+/pABAwYksaqGYUSyYMGC7yX1\nSOSaVIpJNZLknGuS376kh4CHAIqKijR//vyk1M0wjFicc18kek0qR3PWOOd6A4T+rg3IswroE7G/\nXeiYYRhtjFSKyQvAmND2GOD5gDzvAv2cczs657KAU0PXGYbRxkjW0PBk4C2gv3Pua+fcr4GbgVHO\nuWXAyNA+zrltnHMzACRVAJcAs4CPgackLUlGnQzDaF6SYjORdFqcU4cE5P0GGB2xPwOYkYx6GIbR\ncpgHrGEYScHExDCMpGBiYhhGUjAxMQwjKZiYGIaRFExMDMNICiYmhmEkBRMTwzCSgomJYRhJwcTE\nMIykYGJiGEZSMDExDCMpmJgYhpEUTEwMw0gKJiaGYSQFExPDMJKCiYlhGEnBxMQwjKSQUjFxzvV3\nzr0fkX5yzv2+Vp6DnXMbIvL8KZV1MgwjNaR03RxJnwKDoXop0FXA1ICscyUdlcq6GIaRWpqzm3MI\nsEJSwov7GIbR+mlOMTkVmBzn3FDn3IfOuZnOud2bsU6GYSSJZhGT0AJbxwBPB5xeCPSVNAi4G5gW\np4zznHPznXPzv/vuu9RV1jCMRtFcLZMjgIWS1tQ+IeknSZtC2zOATOdc94B8D0kqklTUo0dC6ykb\nKeDNN2HoUCgogH794B//aHxZH38MV1wBv/41PP88VFYmr55G89EsC5cDpxGni+Oc2xpYE1rcfF+8\nwK1rpnoZjWDePBg1CkpK/P7y5XDBBfDDDzB2bGJlPfGEv7a8HCoq4KmnYL/9YOZMyGiuf51GUkh5\ny8Q5lw+MAp6LOHaBc+6C0O6JwGLn3AfA34FTJSnV9TIazzXX1AhJmJIS+POfvSg0lI0b4cILYfNm\nLyQAmzbBW295UTHaFinXfknFwFa1jj0QsX0PcE+q62Ekjw8+CD5eVgZr18K22zasnLlzg1sfxcUw\nZQr86leNr6PR/JgHrJEwO+8cfNw56B5j7QpG8i2SeG3Q/PzG1c1oOUxMjISZMAHy8qKP5eXB734H\n2dn1Xz9jBvTp41seGzfGns/Lg9/+Njl1NZoPExMjYQ47DCZO9N2ZjAwoLITLL4cbb4zNu2kTrFgB\npaV+/7334KSTYNUq3y0Kk57uR4ZycuCyy2DEiOa5FyN5mL28g7NmDfzlL/Dii76LcvnlcOKJvstS\nGwleew2eftq3QKZNg/79fUsiPT06b0WFH9l59FF/zjm47jpYuBC2bIktOy0Nxo+HU06B7bZLya0a\nKcbEpAPz/fcweDCsW+dHYVasgHPOgUWL4Prro/NK3g/kqae8gTQtDR5+GK69FsaNg6++8q2VlSt9\nq2LBAnj88WjhmDABevSAqqrYuuTlwb77mpC0aSS1uTRkyBAZTedPf5JyciQvFTUpJ0daty467+uv\nS/n5wXknTZKys6WMDH8sP19yLjYvSF26SJmZDfvOxrJ0qfTrX0t77imdfrq0aFFyyu1IAPOV4Htp\nLZMOzMsvB3c5srPh/ffh4INh8mS4/35Ytsy3SGojwRlnRI/KBOULU1HhR2p++qmmhZKXBxddBN26\nNel2AF/vX/zCjxRVVvpW1tSp3gnuoIOaXr4RHzPAdmB22MF3V2pTVgbbbAO/+Q2cfz7897/efySI\nsrL4w7tBpKVBZiZsvTV06eK7NXvu6UXguuvgyy8bdSvV/P733ugbdsmvqvIOdRdd1LRyjfpxaoPO\npkVFRZo/f35LV6PNM38+DBsW7c2amQl77w2PPQZDhvj/4VNFWMgibSgZGd72ctZZjSszJ6dm5CgS\n5/zxzMzGldvRcM4tkFSUyDXWMunAFBV50ejWzQ/LZmf7LsL06d47NWhEB2qGg3NyGuZXEr6mNlVV\nscbYigoYMwb+85/Y/CUlvm7Tp8e684fp0iX4eG6uzfVJNSYmHZyTT/bDw++840diXnnFDxH37Bk7\n3BvmZz+DBx/01/XuXf93ZGUl1hUK1yvymunTfZ1OP92nXr38cHZtxo6NdajLzfXdtXjiaCQHExOD\njAzYdVdvxwgzerRveQTx+eeweDF06uQNtAUFwbaXMGVliYcV2LABPv3Ub69Z48WluNgbbn/6ydtF\nTjop1pbzxz/6lk1ODnTu7P8efzzcfHNi328kjomJEUhWlndQCxKUzZvhzjv9SNB++3n/lNqtgaZS\nWgr33OO7QU8+Gb9l83StcFvp6XDffd7vZeZM39qaNMnfj5FarBdpxGX33f3/7kHDx5Lv6tx2G6xe\nnZqARhMn+lZPYWFwaIPyct9KCaJ794ZPOjSSg7VMOjBvvAFHHOEjpZ1xRk23IpKBA4Ovrajwrvdf\nf526yGhbtsBdd3mP2iBDb1aWr7/ROjAx6aBMneon7L30ko+UNmWKH91ZvDg63/XXB3d1ysuDRSQ7\nG+64w9suapOe7icDTphQt40lktLSmhnGkWEJ8vO9AA4e3LByjNRj3ZwOiASXXBI9vFpZ6Q2cV13l\nR07C9O7d8FGQ7GwvUuHWQp8+PiwBeCGZNMlP5APYaSdvKA2ap1O7rsXF8NBDcMIJPtasc3DmmV4M\njdaDiUkHZN06n2oj+UDRkdxxR01IxbqIFJLvv/e2jksu8ecuvRT++U8vJBs2+FbJGWf4c2PG+P2M\njPiu/Vu2eAE54gjr1rRmTEw6IIWF8VsbPXtG78+fX39c10gh+e47GD7cu8lPm+YF5fDDYZddvJAc\ndpgXjpkzowUlSEjAd2d2t5WU2gTNEVB6pXNuUWgd4RgfeOf5u3NueWghrr1TXaeOTnY2nH22d+aK\nJC8Prrwy+tiee9bvOXr55V5ISkq8kCxZArNmwXHHeZGIFJJ58/xcn9GjfdfqjDN8XWqTluZtNRMn\nmudqW6G5DLDDJQ2O4+t/BNAvlM4D7m+mOnVo7rzTB0HKzvYtlbCQ1H6xL7+8fpf5efN8VygvzzuX\nhZk1C446yrdaRo3y+cC3is44w9tRNm/2AZMiCYtHVZXvIv3739HnKyvhhRfg6qv9jOb16xO9eyMl\nJBqzINEErAS613H+QeC0iP1Pgd51lWnxTJLHunU+3semTfHzzJsnDRkSHJ8knE4+WSov9/knTIif\nzznpgQd8vpIS6ZBD6i4XpLw8adkyf82mTb4uBQU15zp1khYuTO1z6mjQiHgmzSEmnwPvAwuA8wLO\nTwcOjNh/BSgKyHceMB+Y37dv31Q8PyOCqirp5Zels87yafbsxATlj38MzvO3v/nzDRUS8MGUxo71\n102YEBzQqX9/X2cjObRWMdk29Lcn8AFwUK3zDRKTyGQtk9Rz0UXRkdWCoqzVTunp0ty5/vr584Pz\nrFrlzz/3XMOEJJwOP9xft8suwedzc6UvvmiZZ9UeaYyYpNxmImlV6O9aYCqwb60sq4A+EfvbhY4Z\nLcQHH/jQBJER0+qKngY1fiQHHuiNrRdfHJzv3HO9Ufb44+HuuxtWn9xcb9gNf08QUsMd4YzUkNLH\n75zLd84VhreBQ4FaPpa8AJwVGtXZD9ggaXUq62XUzUsvBQ8HOxf8Mkc6pEWO2oQpKKjZjhzlueSS\nhgmKczXr6Jx7buwolHM+LIIFo25ZUq3lvYA3QusIvwO8KOmlWmsNzwA+A5YDDwMWYK+ZkHyE+b59\n/TyXQYNgzhzv2xEUkSwrK3hk59JLvZCUlEQLiXNw771+tGX8+Jr8s2b5lkllpReUsL9JPDZvrmkZ\njR0L++9fU8eCAj+h78knG/UIjGSSaL+oNSSzmSSHO+7woyG1bQ/PPuv/Bo3EBNkrjj66pszf/74m\nb3jUJkzkKM9NN/ljZWXS/vvXbzMJ55e8oXXuXOm226TJk6XNm1P/rDoatEYDbCqSiUnTKS+XOncO\nfnGHDpWmTfNG106dfMrKktLS4r/sv/1tTdmXXho9/PuXv9SM8lx/fbSQHH98/UICfkTJaD4aIybm\nW9hB+fHH+C7sH38Mxx4L337rl8OQfIjEiRPjl/fww/7vQw/5+TzguydHH+1DQb7/vrerXHedP1de\n7rtGU6c2rL677tqwfEbLYWLSQena1dtAgiK577KL/1tQ4I2lAEuXehuIFL/MSEGJFBLwKwGCFxQp\nMSHJyIBTT21YXqPlsMG0DkpGhnefDwq+fMMNsfnPPrthkegffhh+/nMfZyQsJGGeesrHTNl774YL\nCXgj7g47NDy/0TJYy6QDc/XVfjLdTTf5kAQ77gi33x4cJ6RXLx8T9uijfYiBSLKzo1s477wT/zs/\n+CD+uawsP0KTlubn+qSleQEaPTqx+4qkqsp323JzfQwVI3XYIlwG4Idp4zmERSL5ANKvvOKHkefN\n82EHgrpLiXL44d5WU1YGu+3mFwhryqJZr77ql8XYuNGLys47+xZRuBtnxKcxi3BZy6SDU1kJ773n\nWwGDB9fvReqcN9wuWuRFJFlCAn7hr7lz/faQITB0aOPF5Kuv4Jhjoj13lyzxAvXFFxbWIBWYzaQD\nM3euX1N4+HD/km23XbTnaiSSX1Zi0CAf4+S+++D11xsuJA0J/VhcXJPeeMMHur76aj+qlCgTJ8Z6\n8Uq+lTJnTuLlGfVjYtJBWbfO2yLWrvULWm3a5JesOPRQ/8LV5vzzvRF20SLfZUi0d5xo/qoq+OYb\nP8y8++5+4a9E+PJL310KKne1TdZICSYmHZQpU4Kjy1dW+qDNK1fWxH795BMfwzXe+r6ppLTUu+Nf\nc01i140cGR3NPkxVFRxwQHLqZkRjYtJBWbPG+4LUpqTEz3/ZYw/o0cMP9daOdNbcVFXB7NmJXXPi\niX70JnKZjvx8f7x//+TWz/CYmHRQhg+Pns0bRvItkuJi3yL4/e/rN1jm5noP1VQaNTt3Tix/VpaP\ntH/ttb6bVFTkZyg//nhKqmdgYtJhOfhgH3ukvjWCS0p8yySeUGRn+5bMvHm+vHiLnccjO7v+EZu8\nPD8zOVEKCnz3aPFiePddOOcci3mSSuzRdlCc80GZ77rL2xCGDInv4fr1196vZNtt/QvaqZP/O24c\nfPaZd3orLPRObdOm1b1I+KGH+iVH8/J8i+Hpp2HGDD9KlJ7uneP69fOtnc6dfZ3OPBMussAUrR5z\nWjMAP/LRs6cPbhRJWpq3Mzz5pLddLFzo8+6zT3CLQvJisGJF7LnttvPeqEHdq9osXeq7V3vs4VcV\nNJqXxjitWcvEAHxr4qabors9aWl+f8KEmv2iorqdyZyDyZN9SyUcES0720dCa6iQgM8/apQJSVvC\n/ACNai680HdlbrzRd22GDvWT/gYMSKycffbxfiH/+If39zjwQO+Nal6n7Rvr5hiGEUOr6+Y45/o4\n515zzn3knFvinBsbkOdg59yG0PKh7zvn/pTKOhmGkRpS3fCsAC6XtDAUpX6Bc26OpI9q5Zsr6agU\n18UwjBSS0paJpNWSFoa2NwIfA9um8jsNw2gZmm00xzm3A7AXEDQvdahz7kPn3Ezn3O5xrj/POTff\nOTf/u+++S2FNDcNoDM0iJs65AuBZ4PeSfqp1eiHQV9Ig4G5gWlAZkh6SVCSpqEePHqmtsGEYCZNy\nMXHOZeKFZJKk52qfl/STpE2h7RlApnOue6rrZRhGckmpAdY554CJwMeS7oiTZ2tgjSQ55/bFC9y6\nVNbLaJ0sXOjn+PTp40M4ml9K2yLVP9cBwJnAIufc+6FjVwN9ASQ9AJwIXOicqwA2A6eqLTq/GI2m\nvNwvF/raa94dPyPDz8uZO9ei0rclUiomkt4A6gzYJ+ke4J5U1sNo3dx5pw/+HBlfpbjYr5Xz9tst\nVy8jMWxujtHiPPJIbKCmqiq/CuCaNS1TJyNxTEyMFideUGrnYoNCG60XExOjxTnllOBYKn36+ImH\nRtvAxMRoca6+2q8mGA5PkJPjQxhMmtSwJTKM1oENvhktTufO3j7y7LN+BGennWDMGB+syWg7mJgY\nrYLsbPjVr3wy2ibWzTEMIymYmBiGkRRMTAzDSAomJoZhJAUTE8MwkoKJiWEYScHExDCMpGBiYhhG\nUjAxMQwjKZiYGIaRFExMDMNICiYmhmEkheaITn+4c+5T59xy59xVAeedc+7vofMfOuf2TnWdDMNI\nPqleazgduBc4AtgNOM05t1utbEcA/ULpPOD+VNbJMIzUkOqWyb7AckmfSSoDpgDH1spzLPCEPG8D\nXZxzvVNcL8MwkkyqxWRb4KuI/a+JXWu4IXlseVDDaOW0GQOsLQ9qGK2bVIvJKqBPxP52oWOJ5jEM\no5WTajF5F+jnnNvROZcFnAq8UCvPC8BZoVGd/YANklanuF6GYSSZVK/oV+GcuwSYBaQDj0pa4py7\nIHT+AWAGMBpYDpQA56SyToZhpIaUB5SWNAMvGJHHHojYFnBxquthGEZqaTMGWMMwWjcmJoZhJAUT\nE8MwkoKJiWEYScHExDCMpGBiYhhGUjAxMQwjKZiYGIaRFExMDMNICiYmhmEkBRMTwzCSgomJYRhJ\nwcTEMIykYGJiGEZSMDExDCMpmJgYhpEUTEwMw0gKJiaGYSSFlIVtdM7dChwNlAErgHMkrQ/ItxLY\nCFQCFZKKUlUnwzBSRypbJnOAPSQNApYC4+rIO1zS4I4kJGWUMYlJnM3ZXMd1rGRlS1fJMJpEylom\nkmZH7L4NnJiq72prFFPMUIayghUUU0wWWdzBHUxlKodyaEtXzzAaRXPZTM4FZsY5J+Bl59wC59x5\n8QpoT8uD/p2/s5SlFFMM+FZKCSWcwRlUUtnCtTOMxtEkMXHOveycWxyQjo3Icw1QAUyKU8yBkgYD\nRwAXO+cOCsrUnpYHncIUtrAl5vhmNrOEJS1QI8NoOk3q5kgaWdd559zZwFHAIaH1cYLKWBX6u9Y5\nNxXYF3i9KfVq7eSRF3i8iipyya3e/5IvmcMcCinkSI4kn/zmqqJhJEzKujnOucOBPwLHSCqJkyff\nOVcY3gYOBRanqk6thYu4KEo0AByO7dmefvQDYDzj+Rk/Yyxj+Q2/oTe9mcvclqhu++b77+G992Dj\nxsZd/803MHs2LF8eP8+GDVAS+Aq0K1JpM7kHKATmOOfed849AOCc28Y5F17hrxfwhnPuA+Ad4EVJ\nL6WwTi3OUpZyC7dQRhngRSSPPLrSlR70YG/25izO4lZupZRSiilmY+hzDMdQSml1WYtYxBVcwcVc\nzKu8ighs/BlBlJbC6adDnz5w8MHQqxdccw0EN6BjqayEc8+FnXeGk0+GQYPg0EOhuLgmz/vvw157\nQY8e0KULHHUUrF2bkttpFUhqc2nIkCFqi5SqVD3VU05ORHwylalc5Vbvpykt6nz400mdNFMzJUl3\n6S7lKlfpShdC+crX6TpdVapq4btsI1x4oZSbK3n58CkzU7r4Yqm8vP7r//rX2OuzsqSzzvLnv/1W\n6tQp+nxGhrTrrlJV6/+NgPlK8L00D9hm5EVeZDObY1oQ5ZSzmc3V+1VUxS1jC1v4lm+5kivZzObq\n0Z9iipnGNF7l1dRUvj1RVgaPPQabN0cfLy+H++7zrY0vv6y7jNtui72+rAwmT/Z/J070fyOpqIBl\ny+CRR5p+D60QE5Nm5Bu+oZzyRl9fTjnDGc5sZpMRYDsvppinebopVewYlJR44QhCglWrfBcoHlVV\n3tYSREWF70J98glsiR2xo6ICLrkEHnzQf9dbb8GECXDXXfDtt4nfSyvCxCQFVFHFLdxCT3qSTjqD\nGMSrvMp+7EdaAo/c4aq3s8nmPu6jM53JJjvqXJg00uKOFBkRdO4MGXUMZFZWwrx5sGSJ354/H045\nBYqK4A9/gJdfhvT04Guzs6GwEPbfH/LjjL6VlcHvfgd9+8KwYV5MrrrKt4hmzAi+pi2QaL+oNaTW\nbjMZp3HKU16UvSNXuXpbb+sYHRN1Lkc5ylNeoJ0k0raSoxw9qkclST/pJ+UrPyZ/rnL1nt5r4btv\nI9S2ZzQ0ZWX5a2vbS8Jp9919+Rs3StttJzmXWPn5+dILL0gff9yij4dG2ExaXBgak1qzmBSrOMqY\nGvkJi0YXddE22kY/0890ra7VL/XLauFwckpTWoyRFqHO6qxyeePgLM1SvvJVqELlK185ytHtur2F\n774NccABjROTcAoSidxc6ZFHar7jm2+kvn0TLzs/35f1i19IP/7YIo+nMWJi3ZwkM4c51cO+tSmh\nhCqqWB/6jGUsi1nMi7xYbZQVooqqwGHeDWxgOd6fYQQjuJu7GcAABjOYp3iKy7gsdTfW3rjhBsjN\nrT9fPBT6fZyDtDTfpTnkEBgzxh9/5hk/LNwYO0hxsTfuzpsH55zT+Do2N4mqT2tIrbVl8ppeU6Yy\nA1slQZ+u6qoc5TQ4P0JZytKf9WcdrIOruzpOTnnK0626taUfQdti+nSpXz/fGigokNLSGtdKycyU\nevWSVq3y5c6dK+XlxebLyEj8O7KypA0bmv3RYN2c5FOhCr2rd7VQC1WpyjrzbqttExKGdKWrkzol\ndA1C2coOFK0c5WiN1jTTk2lHVIZ+1yOPTNzGESkUJ57oyzniiPh5Lr44WGjipZycGpFqRhojJtbN\nqYPXeZ1t2IYRjOAgDqIvfXmXdwPzVlHFKlYlVH4mmVEerQ2llNLAIeZMMvk3/064vA5PWpoftTnq\nKMjMbFwZFRXw//6f316xIjhPXp73mv3VryAnBwoK6v++nj2hd+/G1amZSVk8k7bOd3zHaEZXhwkA\n2MQmRjGKr/maAgqi8gcN1dZHGWWNui4eDkchhUkrr8NQXAwHHQRLl/ph24wMLw6Jkhb6v/mAA7yg\nVNYKJ1FRAf37w403+qHmb76Bd9/1Dm61HeDA1+Phh71dpg1gLZM4TGZyYGyRSip5judijjscW7FV\nQt9RRVWj45cEiVAGGRzCIY0qr0Nz443w0UewaZPfb4yQZGbCCSf47Wuu8QbZSBHIy4OxY32evn3h\n+OPhsstg4MAaY25QmQMHJl6XFsLEJA5rWBMYc2QTm3iKpwJd3h/jMdKJdmZKI41d2IVMEm8+x2u1\npJHGOZxDDjl0Cn26053ZzCaLrIS/p82zZg1cdx2MGOG9S5ctS+z6f/4z2Fu1oRQUwI47ei/WcH3S\n0rxjW3i058IL4fXX4bXXvIfspk2wbh1cemnNCFBtqqrg/vsbX6/mJlEjS2tIzWGAnaM5gY5hYQPo\nlboy8LrZmq09tadylauBGqgX9aIkabmWK0tZCRlaM5QRtw4jNEI/6AdN1VTN1mx9rI/1iT7peBP9\nVqyQunWTsrNrjJz5+X5EpaE01BckM1O67z5pypQap7f0dKlrV+mrr3xZZWXSTjvFXpud7Y2pQf4q\nQ4f60aSg7zzmmNQ8t3rARnOSR6UqNUqj4s7gzVWuilWcUJlH6siEBCVDGXG/P13p+kbf6EN9qH7q\npzzlKV/52l7b6229naKn0gr55S+Dh1t33bXu65YulY49Vios9Ckjo34hmTq15vo336wRlKws6ZVX\nas4tXiz16BF7fWZmcNmDBgV71ObkSDfckJrnVg8mJkmmTGXqoi6BL3OBCrRUS+Neu1AL9Qf9Qf+j\n/9F0Tdcf9AftqB2Vr3wdraM1XdO1j/aJKXeYhmm6pmuYhtUpNJnK1Jt6U93ULeZcoQq1Tuti6rRG\nazRLs7RYi1P52JqXzp2DX9CMDGn9+uBrvvoqdnjWufg+IJFCsm6dtHKl337zTalLF+mxx/z+pk1e\npCTpww9jBSUoZWdL114rHX54dMvFOd/iWbs2pY8vHiYmKeAwHRb4MucqV5u0KfCaS3VpdZwRVOMi\nj9BxOk5lVWWSpPVV67Wv9q3ON0IjVFJVIkkqqSrRCI2IKybZytZETVSBCgLrdo/uqa5Plap0uS5X\njnLUWZ2Vpzzto330nb5L/QNMNdttF/ySZmVJW7YEX7PnnvEFKCurbiHZay9p++2lzz+vOSZ5IfnF\nL6SePaUlS/yxhghKt24+9smcOb5L07mzF5gjj5SWLUvlk6sTE5MU8Lbejplrk6c8/VF/DMx/ha6I\nKwCRQvLJd59IqhGUSCEJn6tLUP6qv+o23RbX43acxlXX6Qk9EWN7yVSmRmlUKh9d8/DXv8a2MrKz\npV/9Kjbvt99Kv/lN/Bc7N1c6/vjoY6efXnP90KE1xyMFJSwk4XO9ekk//ODP3X57/O9zTrrgAm+z\nKSz03aacHGn8+FQ/tXoxMUkRr+gVDdIgpSlN3dVd/6v/VYlKYvKt0IrACXq1heTmN25WxvUZenLx\nk5K8oISF5MH5DyptQpoeWvCQpPiC0lVddYJOCJxUWKACzdCM6nrtpb0C65StbK1VyzSjk0ZFhTRm\njH8JO3f2gjB8uPTTT9H5vv5a6t69ftvIgAHRXrB9+0qrV/syZs+O7or07Sv9/e/SvvtGl/G3v/n8\nW7ZIw4bF/67sbF+n2l63eXnSzJmJP4vJk6U99pC22sq3bD78sNGPtVWJCTAeWAW8H0qj4+Q7HPgU\nWA5c1ZCyW2puzlf6SofpMGUoQ+lK13AN12f6rPr8Hboj8KXNUpY2Vfku0SMLHxHjEeOJEhTJC4kb\n78R45MY7PfvRs5Kkb/VtYLnpSlemMqPm9+QqV0M1NMr1fwftEHh9nvK0TC3XlE4qX34pzZgRf+r+\n+efXLyTa2RNXAAASGElEQVTx0q67xheU2ilSSI48sv6y49lpjjwysfu//fboFppzfoQo3OVKkNYo\nJn+oJ086fh3inYAs4ANgt/rKbgkxKVWp+qhPlC0kTWnqpV7VozrjNT7wpUXoeT0vSVqzaY12u3e3\nKEH539f/V+NeHlctJIxH+zy0j9Zv9gbEx/RY3HIzlKGhGqo9taf20B66VbdqszZH1f1iXRzYHeqp\nnvXON2o37LBD44QkSFCmTQvO86c/+fMNFZK60tChDb+3LVuCh5bT0qSTTmrU42qLYrI/MCtifxww\nrr6yW0JMntEzKlRhzAtZoAJdrat1nI6LG8ck3KUIB4OuLSi1U6SQPKknlaGMuOWGWyin6bS4rYzV\nWq2ttXV1CyZd6cpTXrXAdQiGDGnay52WJj3zjC9rzZrgFsWsWf78okWND74U7ubceWfD7+3TT+P7\nqfTt26jH1RrF5AvgQ+BRoGtAnhOBRyL2zwTuiVPeecB8YH7fRj6gpnCzbo77Ujc07EBtQdn29m1j\nhKT/3f0TEpLIVlInddIKrQis//f6Xtfreg3TMJ2jc/SBPmi2Z9eizJsnXXeddNFFic3WrS0kjz/u\ny9u4MdrYGpm23176LNTtfeutWEGJF52tdp6BA71Rt6Fs2BC/6zVsWKMeW7OLCfAyftGs2ulY/Jo4\n6XiX/RuBRwOub7CYRKaWaJnM0IzAYdhEP6M1urrMo/91dIyYjJk6pvr8/to/obLTla5zdE6zP5tW\nSXm5N0ZGvliNsZnUJyS77BLdSqlLUHbd1Y/a1PVdJ5wglcQa9+vl3HNjxSovzw85N4JW1TKJ+hLY\nAVgccLzNdHMqVKGBGhjlwZqhjLijN0GfERqh4ipvX3lowUNRNpIgo+z6qvWBjm11fXbRLnXexwZt\n0CW6RFtpK3VTN12gC/SDfkj582t2zj47+IXNykpMVE45pabM4cOjz4WNrbNnR7/I229fMzR8yy2x\nglFXzJSxYxt3v6Wlfpg5J8ennj2lSZMa/fhalZgAvSO2LwWmBOTJAD4DdowwwO5eX9ktNZqzXut1\nkS5SV3VVZ3XWKTqlwZHS6hKSTjd1Ut6NeUkRlIN1cNz6V6hCgzRI2cquzp+lLA3QgOrYsu2G8Fyd\noNSvX42NIZ6Le2TLIxyH9ZlnaoSotk0jUlDChtiKiviBkoJSQYEPJi15n5i//c2XNXduwxfu2rzZ\nG4orm2ZYb21i8g9gUchm8kJYXIBtgBkR+UYDS0OjOtc0pOzWFLbxAB1Qr80kS1n6scr/g/znB/+M\nGbX5oeQHrd64OmaUZ/qn0yVJX+tr5Su/QbaZnbVzXGGI11UrVKGe03PN9syahfT0+C/tf/7jX9oJ\nE6QnnpCuuqrucIpDhkQLyl13+e0tW6RLLokeNr7+er9dUSGddlrd4hFp58jNlQ480F83c6bvouTk\n+FZMfr6P4tZEgUiEViUmqUwtISZVqgocRl2jNfq5fh71v33Q51W9Kkla+eNK7XjnjjGjNlL0KM9u\n9+6mNZt8CMaP9JHO1/lxJ/1FfnKUo2t0jbYo1pX8L/pL1NB25OfP+nNqHlxL0b9/8Auclha8/OeX\nX0p/+Uv8Fs3IkdGR4rdskUaP9ucGDKgRFKlhQgI+Qn5YMDIypKIi6ZNPgkeC8vOlJ5+MrXeKMDFJ\nAT/qR52pM5WlLKUpTSM0Qp/q06g8S7REAzSgzpe8UIV6Q29I8oJy0lMnVQvJU3oqapTnl0/+slpI\nXtJL9QpV7U+mMtVFXfS6Xo+q55N6Mu7w9hN6ItWPsnl5993g1ka45RDEJ58Ez/XJzfUve7iFEikk\n4RQWlIYKCcTabtLSvEdsPCNtoo5sTcDEJMlUqUp7a+8oo6uTUzd1qzZabtAGdVXXBhliIwUlTHj4\nN3LYOMxMzUxYSGp/X+RkxFKVahttE+N411M9A6cHtHmWLPGjL506eTtJZAiB2tx1V/yuzu6715zr\n3l3q0yc2T3q6n0G8227e9rHVVnV7yobFo/axvLz4Q8hHH91sj87EJMm8oTcCbQx5yqte8OpBPRiz\nel9DBaW2H0mkoDRVSMLf9aSim8Zf6AuN0AhlhD7DNCyub0qHYfny2NnCkal797pFITxh76OP/PIZ\n8+Z5g+lVV9VtuwlKubnB/jD5+dJzzWfXaoyYWEDpOviUTxGKOV5CCR/yIQCf8zkllDS4zI1s5BAO\nYQhDeJu3o8I/llLK0RzNvuzLO7xDBQ2LRZpGWmAYySqq2MjGqGN96csrvMJmNiNkaxMD/OMfPpB0\nPEpLfejFqthnDNS88rvu6hPAF1/AnXfGBpUOk5XlQzqW1lqdwDl/3aWX+v2yMh8L9qST4LjjEruv\nZsbEpA52Z/fA43nkMYQhAOzDPhRQwCY2NbjcUkp5kzfJIYdKKqOWraiggjd5s8FlFVDAhVzIPdzD\nZqIjnFdSyShGBV6XSxNWs2tvzJtX9/mSkvhCEubBB/2C6H/9q99//nkvMEFkZMAuu3ihWbmyRlDy\n8nwc29/+Fk48EZ59Fn78EUaO9KsDtnYSbcq0htScNpP9tX9UdyNd6eqhHlovbzwtV7l21+4JOa8h\nlK983a7bAw2iiXx6qqdKVapjdEzUCn/5ytd1uq5ZnlObp7YxtbHJuZpRnbFj4+cbPdq7y69fL115\npZ+E2L+/dNttPoZsKwCzmSSfjdqoi3SRClWobGXrGB2jlVoZlWeDNuhQHdpgAchQhsZpXODEwERE\nqad6Vs+xqVSlntNzOkknaYzGxIzkGHVwyy2J2zbipYd8HJrqZUeDjK6ffNKy99sAGiMmzl/Xtigq\nKtL8+fNbuhpRlFBCf/rzNV/Xm3cXduFzPm/0mjlZZHEmZ/IgD8YsrWE0gh9+gJ/9zC890VSeeMJ3\nUQoLg+0lmZl122daCc65BZKKErnG1s1pIkL8wA9kksl85rMbu9V7zXKWN1pIutCFN3mTR3jEhCRZ\ndOsG77wDRx7p7Rnp6T7FI94Ke+npcMwxXjDiLfvZtSu88YZfRbCdYWLSBCYzmW3Yht70pgtduJVb\neZ/32Y7tUvadG9nI9myfsvI7LDvtBNOnQ3m5X9HvppuC86Wnw//8D+y8c/TxtDT417+8ETYjwy+s\nlZMTncc53wo66igvKr16+ZX/9tnHL87V1km0X9QaUmuYmzNLs2L8S/KUp9/pd1qt1SpSUYNtH4l8\nnJw2aENL337754orgmf3ZmXVhBh49VXpvPO8V22kq73kwwgce2xNbNq6ZgqH/Utee63ZbzMemAG2\n+ThABwS+7JGLc12tq+POhWnsp5d6tfCddxBqhxsIp86dvWNaQ7n7bh+tviHG2/32S939JEhjxMS6\nOY3kMz4LPJ5GGt/xHQA3ciPv8R6FFEatG5xFFplkkk02AOmkk0ceB3Mw2WSTFudncTge47Ek34kR\nyODBwXaPsjJvrG0IDz8MV17p1x5uCEuWNLx+rRATk0ayN3sHLiyeTjq96V29P5CBfMRHnMIpFFJI\nD3pwGZexjGVcwRUMYxi/5tcsYAGv8RpLWMIzPMPrvM5RHEUmmWSRRSc6MYUpHMERzXmbHZexYyE7\nO/pYTg4MGwb9+tV/fVUVXH21d3hrKDvskFAVWx2JNmVaQ2oN3Zz39F6gzeQO3ZHU71mndVqhFapQ\nRVLLNRrAggV+TRznvE3jgguk4gauL71+ff2Bl2pP8Js2LbX3kwBYN6f5GMxg/sN/OIRD6ExnBjCA\nh3mYS7k0qd/TjW7sxE42DNwS7L23d7UvL/dDufff713eI3nvPTj4YN9q2XpruPlm719SUAC5caYs\n5OTAqFHQqZMf+end27vjH3tsym8pldjcnCZQRBEv83JLV8NINfF8TpYuhYMOgk2heVlr1sANN8BX\nX8G993p7yY03Rnd18vJg0iQ/aa+qCrZs8aITz3elDWEtk3ZABRU8zdOcyZlcxmUsoW0b8toMN98M\nm6MnV1JSAhMnem/acePguuugS5eaFsgDD9TM/k1L8+LSDoQESJ07vXPuSaB/aLcLsF7S4IB8K4GN\nQCVQoQa48LZGd/qWoowyRjKS93iPTWwinXSyyOJBHuRMzmzp6rVvBg2CRYtij3fuDLNmwc9/7vfb\nYAukMe70KevmSDolvO2cux3YUEf24ZK+T1Vd2jrf8i2TmMRa1nIIhzCSkdXDx//iXyxkIcV49+xK\nKtnMZi7gAk7gBPLJb8mqt28GDvTDubXDE2zZAjvuWLMfboG0c1LezXHOOeBkYHKqv6s98gqvsDM7\ncy3Xcgu38Et+yaEcWh0DZQpTqoUkkgwy+C//be7qdiyuuirWZT43F04+GXr2bJk6tSDNYTP5BbBG\n0rI45wW87Jxb4Jw7rxnq02aooIJTOIUSStjCFgA2sYm3eIvHeRyAQgoDrxWyVkmqGTgQZs6EPfbw\n3Zf8fLjoInjkkZauWYvQpG6Oc+5lYOuAU9dIej60fRp1t0oOlLTKOdcTmOOc+0TS6wHfdR5+vWH6\n9u3blGq3GRaykDJip6uXUMITPMFv+S3ncz4zmBETOjKffPZjv+aqasfloIO83aS83BtZ24hNJBU0\nSUwkjazrvHMuAzgBQjEOg8tYFfq71jk3FdgXiBETSQ8BD4E3wDah2m2GdNIRwbeaEfrpRjKSy7iM\n27iNDDJwOLLIYiYzzTelOYkXcqADkepuzkjgE0mBEYOcc/nOucLwNnAofuFzA9iLvehEp5jj+eTz\nG35TvX8DN7Cc5TzAA0xhCqtZzWBiBs4MI6Wk2mntVGp1cZxz2wCPSBoN9AKmehstGcC/JL2U4jq1\nGdJI43meZyQjqwNPp5POsRzLaZwWlXdbtuV0Tm+hmhpGisVE0tkBx77Bry+MpM+APVNZh7ZOEUWs\nYhXTmMb3fM8whlmrw2iVmDt9GyCffGt1GK0ec6c3DCMpmJgYhpEUTEwMw0gKJiaGYSQFExPDMJKC\niYlhGEnBxMQwjKRgYmIYRlIwMTEMIymYmBiGkRRMTAzDSAomJoZhJAUTE8MwkoKJiWEYScHExDCM\npGBiYhhGUjAxMQwjKZiYGIaRFJokJs65k5xzS5xzVc65olrnxjnnljvnPnXOHRbn+m7OuTnOuWWh\nv12bUh/DMFqOprZMFuPXxYla58Y5txs+Mv3uwOHAfc65oEVcrgJekdQPeCW0bxhGG6RJYiLpY0mf\nBpw6FpgiqVTS58By/OJaQfn+L7T9f8BxTamPYRgtR6qi028LvB2x/3XoWG16SVod2v4Wv45OIJHL\ngwKlzrn2uFhXd+D7lq5Eimiv99Ze76t/ohfUKyYNXE+4yUiScy7usp+Ry4M65+ZLKoqXt63SXu8L\n2u+9tef7SvSaesWkvvWE47AK6BOxv13oWG3WOOd6S1rtnOsNrG3EdxmG0QpI1dDwC8Cpzrls59yO\nQD/gnTj5xoS2xwBJa+kYhtG8NHVo+Hjn3NfA/sCLzrlZAJKWAE8BHwEvARdLqgxd80jEMPLNwCjn\n3DL8Iuc3N/CrH2pKvVsx7fW+oP3em91XCCfFNVMYhmE0GPOANQwjKZiYGIaRFNqMmDTVdb+t4Jwb\n75xb5Zx7P5RGt3SdmoJz7vDQ77LcOdeuPJydcyudc4tCv1PCQ6mtBefco865tZG+W42Z6tJmxISm\nu+63Jf4maXAozWjpyjSW0O9wL3AEsBtwWuj3ak8MD/1ObdnX5HH8uxNJwlNd2oyYJMF132h+9gWW\nS/pMUhkwBf97Ga0ISa8DP9Q6nPBUlzYjJnWwLfBVxH481/22xO+ccx+Gmp9teSZ1e/xtIhHwsnNu\nQWi6R3uiwVNdwqRqbk6jaC7X/ZamrvsE7gduwP9DvQG4HTi3+WpnJMCBklY553oCc5xzn4T+l29X\n1DfVJUyrEpMUu+63Ghp6n865h4HpKa5OKmlzv00iSFoV+rvWOTcV361rL2KS8FSX9tDNaajrfpsg\n9MOFOR5veG6rvAv0c87t6JzLwhvKX2jhOiUF51y+c64wvA0cStv+rWqT8FSXVtUyqQvn3PHA3UAP\nvOv++5IOk7TEORd23a8gwnW/jXKLc24wvpuzEji/ZavTeCRVOOcuAWYB6cCjoakW7YFewFTnHPj3\n6F+SXmrZKjUO59xk4GCge2h6zJ/xU1uecs79GvgCOLnecsyd3jCMZNAeujmGYbQCTEwMw0gKJiaG\nYSQFExPDMJKCiYlhGEnBxMQwjKRgYmIYRlL4/3N7IIdpK1d4AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f7a2166b710>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "gmm = MultivariateGaussianMixture(n_components=3)\n",
    "gmm.fit(x_train)\n",
    "p = gmm.classify_proba(x_train)\n",
    "\n",
    "plt.scatter(x_train[:, 0], x_train[:, 1], c=p)\n",
    "plt.scatter(gmm.mu[:, 0], gmm.mu[:, 1], s=200, marker='X', lw=2, c=['red', 'green', 'blue'], edgecolor=\"white\")\n",
    "plt.xlim(-10, 10)\n",
    "plt.ylim(-10, 10)\n",
    "plt.gca().set_aspect(\"equal\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
